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Abstract 

We use leading-order anisotropic hydrodynamics to study an azimuthally-symmetric boost- 
invariant quark-gluon plasma. We impose a realistic lattice-based equation of state and perform 
self-consistent anisotropic freeze-out to hadronic degrees of freedom. We then compare our results 
for the full spatiotemporal evolution of the quark-gluon plasma and its subsequent freeze-out to 
results obtained using U-ld Israel-Stewart second-order viscous hydrodynamics. We find that 
for small shear viscosities, LKyjs ~ 1, the two methods agree well for nucleus-nucleus collisions, 
however, for large shear viscosity to entropy density ratios or proton-nucleus collisions we find im¬ 
portant corrections to the Israel-Stewart results for the final particle spectra and the total number 
of charged particles. Finally, we demonstrate that the total number of charged particles produced 
is a monotonically increasing function of 47rr//s in Israel-Stewart viscous hydrodynamics whereas 
in anisotropic hydrodynamics it has a maximum at A-Kyjs ~ 10. For all A'Kyjs > 0, we find 
that for Pb-Pb collisions Israel-Stewart viscous hydrodynamics predicts more dissipative particle 
production than anisotropic hydrodynamics. 


1 



I. INTRODUCTION 


In recent years there have been signihcant advances in our understanding of the theory 
of relativistic hydrodynamics and its application to describing the spacetime evolution of 
the quark-gluon plasma created in ultrarelativistic heavy-ion collisions. Early works applied 
relativistic ideal hydrodynamics [1-3] and later works developed and applied second-order 
viscous hydrodynamics [4-35]. In recent years, the framework of anisotropic hydrodynamics 
was developed in order to better account for the large momentum-space anisotropies gen¬ 
erated in relativistic heavy-ion collisions [36-53]. This framework has been shown to more 
accurately describe the evolution of systems subject to boost-invariant and transversely 
homogeneous 0-l-ld flow than traditional viscous hydrodynamics approaches [47, 49, 51, 53- 
56] and has recently been shown to best reproduce exact solutions of Boltzmann equation 
subject to H-ld Gubser flow [57-59]. In its latest form, leading-order anisotropic hydro¬ 
dynamics allows for multiple local momentum-space anisotropies in the argument of the 
non-equilibrium distribution function [42, 48, 59]. In addition, it has been shown that it is 
possible to account for non-spheroidal/-elhpsoidaI corrections via a modihed shear correc¬ 
tion resulting in “viscous anisotropic hydrodynamics” [47, 53]. With these improvements, 
anisotropic hydrodyanmics has been shown to work extremely well when compared to exact 
solutions of the massless and massive Boltzmann equation in relaxation time approximation 
[47, 53-56]. 

Despite this promise, turning anisotropic hydrodynamics into a practical phenomeno¬ 
logical tool for use in modeling heavy-ion collisions requires two additional fundamental 
ingredients to be implemented: (1) a realistic lattice-based equation of state (EoS) and (2) 
self-consistent anisotropic freeze-out to hadronic degrees of freedom. For the EoS, it is not 
obvious a priori how one enforces thermodynamic relations in an anisotropic system. As we 
will show, it is possible to impose the EoS as a relation between the isotropic energy den¬ 
sity and pressure. For anisotropic freeze-out, we determine the freeze-out hypersurface by 
specifying a critical energy density and then we use the leading-order anisotropic distribu¬ 
tion function to compute particle spectra using “anisotropic Cooper-Frye freeze-out”. This 
method includes the leading dissipative corrections at freeze-out in a way that guarantees a 
positive dehnite one-particle distribution function at all momenta, thus avoiding the problem 
of regions of phase space where the distribution function is negative. We here restrict our 
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anisotropic hydrodynamic analysis to the inclnsion at leading order of the most important 
local momentnm anisotropies [42, 48, 59]. The generalization to “viscons anisotropic hydro¬ 
dynamics [47, 53] which also inclndes the smaller residnal deviations from the leading-order 
distribntion will be pursned in a fnture work. 

As preparation for implementing the necessary ingredients in a 3-|-ld anisotropic hydro¬ 
dynamics code, in this paper we perform the somewhat simpler task of implementing l-|-ld 
anisotropic hydrodynamics for an azimnthally-symmetric and boost-invariant system. We 
will compare to predictions of a 1-l-ld Israel-Stewart viscous hydrodynamics code in which 
we have implemented exactly the same initial conditions, freeze-out hypersurface algorithm, 
etc. Since the H-ld task is more straightforward computationally, this allows us to virtually 
eliminate systematic computational errors. We compare our results for the full spatiotem- 
poral evolution of the quark-gluon plasma and its subsequent freeze-out to results obtained 
with l-|-ld Israel-Stewart second-order viscous hydrodynamics. We hnd that for small shear 
viscosities, Airri/s ~ 1, the two methods agree well for nucleus-nucleus collisions, however, 
for large shear viscosity to entropy density ratios or proton-nucleus collisions we hnd im¬ 
portant corrections to the Israel-Stewart results for the hnal particle spectra and the total 
number of charged particles. We demonstrate that the total number of charged particles 
is a monotonically increasing function of Airri/s in Israel-Stewart viscous hydrodynamics 
whereas in anisotropic hydrodynamics it peaks at Airri/s ~ 10. For all Airri/s > 0 we hnd 
that for Pb-Pb collisions Israel-Stewart viscous hydrodynamics predicts harder spectra and 
more dissipative particle production than anisotropic hydrodynamics. 

The structure of paper is as follows. In Sec. II, we specify the conventions used in the 
body of the text. In Secs. Ill and IV, we present the setup and details for the conformal 
anisotropic hydrodynamics dynamical equations. In Sec. V, we discuss how we implement 
a realistic EoS in the context of anisotropic hydrodynamics. In Sec. VI we discuss how 
to implement anisotropic Cooper-Frye freeze-out in the context of leading-order anisotropic 
hydrodynamics. In Sec. VII we present our numerical results for Pb-Pb and p-Pb collisions 
at LHC energies. In Sec. VIII we summarize our hndings, state our conclusions, and present 
an outlook for the future. Finally, in six appendices we present details behind many of the 
derivations and results presented in the body of the text. 
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II. CONVENTIONS 


In this paper, the metric is taken to be “mostly minus” such that, in Minkowski space 
with = (t, x, y, z), the line element is 

ds'^ = g^udx^dx'' = dt^ — dx^ — dy^ — dz^ . ( 1 ) 

Since we deal with a system that is boost invariant along the beam-line we transform to a 
new variables dehned by r = \/t^ — z"^ as the longitudinal proper time, and ^ {z/1) 

as the longitudinal spacetime rapidity. Also, since the system is cylindrically symmetric 
with respect to the beam-line it is convenient to transform to the polar coordinates in 
the transverse plane with r = \Jx‘^+y‘^ and (j) = tan“^(j//x). The new set of coordinates 
x^ = (r, r, (j), g) dehnes the polar Milne coordinates. 

We also mention that the notations | and | 

denote symmetrization and antisymmetrization, respectively, and = A^^A"^ where 

= aL^A^^ — A^^Aa/ 3/3 with . The four-index projector A^^ projects 

out components of a rank-2 tensor which are traceless and transverse to the flow velocity u. 

III. ANISOTROPIC HYDRODYNAMICS SETUP AND BASIS VECTORS 

For the purposes of this paper, we assume that the system is boost invariant and az- 
imuthally symmetric so that we can apply H-ld hydrodynamical evolution. In the case of 
anisotropic hydrodynamics we use an ellipsoidal form for the local-rest-frame (LRF) one- 
particle distribution function. This form can be obtained by introducing an anisotropy 
tensor of the form [42, 48, 59] 

+ e'" - , ( 2 ) 

where is the four-velocity, is a symmetric traceless anisotropy tensor, <F is a parameter 
associated with the bulk viscous correction. In the LRF, the anisotropy tensor is = 
diag( 0 ,^ 2 ,,^ 2 ) with ^^lrf = = 0 [48, 59].^ As mentioned above, the held 

<h accounts for bulk viscous effects. In the case of a massless (conformal) system, one can 

^ The labels x, y, and z do not imply the cartesian components, but are merely labels for spacelike directions. 
In this work, the vectors corresponding to x, y, and z will be mapped to the radial, azimuthal, and rapidity 
directions, respectively. 
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take <k = 0. We will assume from here on that, even in the non-conformal case, one has 
<h = 0. This assumption is analogous to assuming that the bulk correction to the pressure 
can be neglected relative to the shear correction in second-order viscous hydrodynamics. 
Bulk viscous effects will be included in a future work. 

In general, can have off-diagonal elements; however, based on the proper-time evo¬ 
lution of the shear-stress tensor components in simulations, one can expect that the off- 
diagonal elements are much smaller than the diagonal ones and can therefore be considered 
perturbatively [60]. Importantly, however, for l-|-ld expansion, similar to viscous hydrody¬ 
namics where it suffices to include only the diagonal contributions to the shear tensor [nl., 
TT^, and TTp, in anisotropic hydrodynamics one only needs the diagonal components of the 
anisotropy tensor to describe a 1-l-ld system. For a discussion of 1-l-ld second-order viscous 
hydrodynamics, see App. D. 

In a general frame one can expand the (diagonal) anisotropy tensor in covariant form 


+ iyY^Y'' + , (3) 

with 

e, = 0. (4) 

The orthogonal basis vectors and Z^ reduce in the LRF to unit vectors in the 

t, r, 0, and 2 ; directions. The basis vectors are explicitly dehned in App. A. They obey the 

normalization conditions = 1 and X^X^ = Y^Y^^ = Z^Z^, = — 1, therefore, = 1. 
In lab frame, and are unit vectors which point in a mixture of the r and r directions, 
Y^ is a unit vector pointing in the azimuthal direction, and Z^ is a unit vector pointing in 
the spatial rapidity direction. 

Using the tensor one can construct an anisotropic distribution function 

, (5) 

where feq{x) = l/[exp(a;) -|- a] with a = —or 0 for Bose, Fermi, and Boltzmann 
statistics, respectively. Above, A is a temperature-like scale that can be identihed with the 
temperature, T, only when = 0.^ 

^ We assume herein that the chemical potential is zero. 


f{x,p) = f, 


eq 
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IV. 1+lD ANISOTROPIC HYDRODYNAMICS 


We will obtain the anisotropic hydrodynamics equations by taking moments of Boltzmann 
equation in the relaxation time approximation (RTA). The Boltzmann equation in RTA is 

W = ^(/-/eq), (6) 

^eq 


where is the covariant derivative which becomes the ordinary derivative in flat space- 
times and Teq is the relaxation time. The right-hand side of Boltzmann equation is the 
collisional kernel which contains all interactions involved in the dynamics. If the particles 
comprising the fluid are massless, conformal invariance requires that Teq is inversely propor¬ 
tional to the temperature, i.e. Teq oc 1/T. In conformal RTA one has Teq = St^/T, where 
fj = 7]/s with r] and s being shear viscosity and entropy density, respectively. 

The general equations governing conformal 1-l-ld anisotropic hydrodynamics for a 
azimuthally-symmetric and boost-invariant system were obtained originally by Tinti et 
al [48]. The resulting dynamical equations give the evolution of six helds ^y, ^z, A, T, 
and 9± where 9± = tanh“^(M'’/M'^) is the transverse rapidity. Using the constraint (4) one 
is left with hve independent parameters. Taking four equations from the hrst and second 
moments of Boltzmann equation one has [48, 59] 


DuE -|- £9u -I- PxDx9± -|- P, 


y 


DxPx 


Px9x + e:Du9j 


sinh 9 i 
r 

cosh 9 I 


+ P. 


~ ^ 2-^ 


Oii 


Pu^j 

a, 


— O'i + 


2n 




eq 


y 

al 


Pz 


cosh 9\ 

T 

sinh dj 

T 
1 


C^x^y^z 


= 0, (7) 

= 0 , ( 8 ) 

= 0, ie{x,y}. (9) 


Above Oj = l/yr+~5, £ is the energy density, and Pi are the spacelike diagonal components 
of the energy momentum tensor (pressures) with [59] ^ 


£ £ gqi^X'jP.i^CX X y Oiz') y (70) 

Px = Peq(X)?{Tx(o:xy C^z) y (H) 

Py = Peq{,X)l-LTy{,Olxy a^) , (12) 

Pz = Peq{X)T-iL{axy a^) y (13) 


^ For non-conformal systems, one has instead Oi = l/vT+”U~+~T. 
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where Seq and Peq are the isotropic equilibrium energy density and pressure, respectively, 
and the special functions 71, Htx, Tiry, and Hi are dehned in App. C. The derivatives 


Ou = dy,U>^ , 

( 14 ) 

Ox = d,X^, 

( 15 ) 

Du = , 

( 16 ) 

Dx=X^Dy,, 

( 17 ) 


are dehned in the 1+ld case in App. B and cxj are the diagonal projections of the velocity 
stress tensor 


a. = = 


cosh 6 I sinh 6 1 


+ 


T 


2^ 

3 


ay = Y^a^^Y^ = j 


6„ sinh 6 i 


( 18 ) 

( 19 ) 


The hfth equation necessary is obtained by requiring energy conservation. In the confor¬ 
mal case, this constraint implies that the hrst moment of the collisional kernel must vanish 
and results in the dynamical Landau-Matching condition 


T = XR'/Ha,,a,). 


( 20 ) 


In the next section, we will describe how to extend this to the non-conformal case. Finally, 
we note that the ads obey the following constraint^ 


1 1 1 

_I_I_= Q 

- 2 ' ^2 ' ^2 ’ 


a 


a: 


at 


which can be used to determine ay as a function of ax and a^ 


( 21 ) 


ay[a 






y/3 


a^al 


at 


at 


( 22 ) 


V. LATTICE-BASED EQUATION OF STATE 


The dynamical equations presented in the previous section were obtained [48, 59] in the 
conformal (massless) limit. In this case, there is no fundamental scale and one has 


Peq t Pideal 5 
^eq t £ideal t 


For a non-conformal system, the right hand side of Eq. (21) would instead be 3(1 -I- $). 


( 23 ) 
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(a) (b) 




FIG. 1. The left panel (a) shows the equilibrium energy density and pressure obtained from 
Eq. (26), scaled by their ideal Stefan-Boltzmann limit. The right panel (b) shows the speed of 
sound (c^ = dPea^/dseq) as a function of temperature. 


with 


^ideal STjdeal 


( 24 ) 


For the quark-gluon plasma the ideal (Stefan-Boltzmann) limit of the pressure is 


TT 


= PSB = ( N. 


45' 


1 + In.Nj 

4 


( 25 ) 


where and Nf are number of colors and quark flavors, respectively. 

In practice, however, interactions induce corrections to the ideal EoS and also result in 
the running of the strong coupling constant which breaks conformal invariance. At low 
temperatures, the quark-gluon plasma behaves non-conformally and one must use lattice 
QCD simulations to determine the EoS. Herein, we employ an analytic parameterization of 
lattice data for the QCD interaction measure, I(T) = Seq — SP^q, taken from the Wuppertal- 
Budapest collaboration [61] 


HT) 


+ 


/o[tanh(/it -F /a) l] 


hi ha 


(26) 


1 + aP l + git + g2p 

with t = T/(0.2 GeV). For Uf = 2 + 1 {2 light and one heavy quarks) the parameters are 
ho = 0.1396, hi = -0.1800, ha = 0.0350, /o = 2.76, fi = 6.79, /a = -5.29, gi = -0.47, 
g 2 = 1.04, and a = 0.01.^ 


® In the original parameterization presented in Ref. [61] the authors used a = 0, however, choosing a = 0 


















The pressure can be obtained from an integral of the interaction measure 

Peq(T) _ dTIjT) 

Ti ^4 

where we have assumed P{T = 0) = 0. Having Peq(T), one can obtain the energy density 
£eq via 

^eq = ^Peq{T) + I{T) . (28) 



Note also, that one can construct the inverse function T{e) straightforwardly. 

In the limit T —)■ oo, the system tends to the ideal limit (24) as expected. The tem¬ 
perature dependence of the resulting isotropic energy density, pressure, and speed of sound 
squared (c^ = dPeq/dSeq) are shown in the two panels of Fig. 1. In what follows, the tem¬ 
perature dependence of the lattice-parameterized £eq and Peq are used when computing the 
anisotropic energy density and pressures via Eqs. (10)-(13) and the inverse function T{e) 
is used where T appears in Eqs. (7)-(9). In addition to using the self-consistent effective 
temperature and lattice-based EoS for £eq and Pgq, one must also specify how the relax¬ 
ation time is determined. In the relaxation time approximation with a general EoS, one has 
Teq = 5?7/4Peq. In the conformal limit, one can use 4Peq = Sgq + Peq = Ts to rewrite this 
as Teq,conformal = 5f]/T. For a general EoS, if one works with fixed fj = r^/s, one has instead 
Teq = bf] {IP £eq/F’eq)/4T. 

Before proceeding, we note that there is somewhat of a inconsistency in our prescription 
for implementing the EoS in anisotropic hydrodynamics. This stems from the fact that the 
factorization (10)-(13), which occurs in the conformal case, no longer holds in the case of a 
non-conformal (massive) gas [42, 49, 56]. However, as we demonstrate in App. E, one hnds 
that for 0.1 ^ Pl/Pt ^ 10 and masses m/T ^ 1 the factorization of the thermodynamic 
variables is accurate to ^ 5%, with the largest corrections being seen for a strongly pro¬ 
late plasma {P^ S> Pr)- Although we have no precise knowledge of the effective degrees of 
freedom and their masses in a QGP in the temperature range considered here, this raises 
some hope that the conformal factorization approach used herein is a reasonable approxi¬ 
mation. However, even with this understanding, there is another complication if one breaks 
the conformal symmetry since, in this case, one cannot naively assume that the parame¬ 
ter <F appearing in the general LO ansatz for the one-particle distribution function (2) is 
gives the wrong high temperature limit. We introduce a small a > 0 which does not affect the parame¬ 


terization in the vicinity of the phase transition, where the original fitting was performed, but guarantees 


that the pressure approaches the Stefan-Boltzmann limit in the high-temperature limit. 
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zero [42, 48, 59]. In a forthcoming paper, we study an alternative prescription for imposing 
the EoS which employs a quasiparticle model with a temperature-dependent quasiparticle 
mass that is tuned to reproduce the lattice EoS [62] as suggested recently in Ref. [63]. Our 
preliminary Endings suggest that this alternative approach agrees quite well with the imple¬ 
mentation chosen here for the evolution of the efiective temperature and shear corrections, 
but that there may be important difierences in the evolution of the bulk pressure correction, 
which, in the current approach, is related to the diEerence between and We leave the 
comparison of these two methods to a future paper [62]. 


VI. ANISOTROPIC FREEZE-OUT 

We now turn to the question of hadronic freeze-out. Our technique will be to perform 
“anisotropic Cooper-Frye freeze-out” using Eq. (5) as the form for the one-particle distri¬ 
bution function. This is difierent than the typical freeze-out prescription used in viscous 
hydrodynamics in which one takes into account the dissipative correction to the equilibrium 
distribution function only at linear order. One immediate benefit of performing anisotropic 
freeze-out using Eq. (5) is that, with this form, one is guaranteed that the one-particle 
distribution function is positive-definite at all space-time points in the plasma. 

In practice, we start from the standard freeze-out integral 

= (29) 

In the integral above, S is the three-dimensional freeze-out hypersurface defining the bound¬ 
ary of the four-dimensional volume occupied by the fluid, is the surface normal vector, 
and is the particle four-current. Due to the presence of momentum-space anisotropies, 
one cannot simply use the momentum scale A when defining the freeze-out hypersurface S. 
Instead, one should use the energy density, from which one can obtain the effective freeze-out 
temperature Tpo = Tes = T{e) using the realistic EoS described in the previous section. 

After identifying S, we use the following parametrization of the freeze-out hypersur- 
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face [15, 16, 64, 65],® 

t = {tq + d{C, 0, 9) sin 9 sin () cosh 
X = d{(, (f), 9) sin 9 cos C cos 0, 
y = (i(C, 0, 9) sin 9 cos C sin 0 , 
z = {tq + d{(, 0, 9) sin 9 sin (■) sinh 


d{(, 0, 9) cos 9 
A 


d{(, 0, 9) cos 9 
A 


(30) 


Based on this one hnds the following expressions for the longitudinal proper time r, longi¬ 
tudinal rapidity <;■, and the transverse (r and 0) coordinates of S 

T = To + d{(, 0, 9) sin 9 sin ( , 
r = d{(, 0, 9) sin 9 cos ( , 


= 


d{(, 0, 9) cos 9 
A 


( 31 ) 


The function d{(, 0, 9) is by construction the distance of the points on freeze-out hyper¬ 
surface to the point (tq, 0, 0, 0) in (r, r, 0, A^)-system of coordinates (tq being initial proper 
time) as can be seen from Eq. (31). The freeze-out surface variables 9 and ( are the polar 
and azimuthal coordinates in this coordinate system. The length scale A is introduced for 
dimensional reasons and hnal results are independent of this quantity. The normal vector to 
the hypersurface is constructed in the usual manner by taking derivatives of the orthogonal 
basis coordinates (30) with respect to the relevant parameters (, 0, and 9 

( 32 ) 


d^T.^ = e^^f,^ — ——dCd^d9, 


d( 00 89 

where e^ajB'y is the four-dimensional Levi-Civita symbol. 

To proceed we use the kinetic dehnition of 

f = j dxp^f{x,t). (33) 

For energy densities below the QCD phase transition temperature (energy density), it is 
appropriate to describe the system as a gas of hadrons; therefore, dx translates to 


d^p 


dx = + l){2gi + l)-^5(p>^ - m^i)2Q{p°) , 

, (27r)^ 


(34) 


® Different parametrizations have been used in the literature [3, 6, 66, 67]. The parametrization used here 
has the advantage that the function d{(^, cj), 6) is single-valued for most (but not all) freeze-out surfaces, 
including the typical one shown in Fig. 8 for which other parameterizations, e.g. using the freeze-out 
proper time TFo(r), are multivalued. 
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where Si and Qi are the spin and isospin degeneracies of the hadron and m* is the hadron 
mass. Putting everything together one has 




P 


with 


P 


,dN 


Mi 


P 




fi{x,p)jfd E„, 


( 36 ) 


( 36 ) 


d^p /. (27r)^ j 

where Mi = (2sj + l){2gi + 1) is the degeneracy factor and /* is the distribution function for 
the particle species i taking into account the appropriate quantum statistics. 
Parameterizing the particle momentum in the lab frame as 


= (m_L cosh y, p± cos ip, p± sin ip, m± sinh y ), 


( 37 ) 


one hnds 


p ■ u = m± cosh(6*_L) cosh(|/ — — p± sinh(6*j_) cos(0 — tp ), 

p ■ X = m± sinh(6'_i_) cosh(j/ — ?) — cosh(6'_i_) cos(0 — ip), 

p -Y = p± sin(0 — ip ), 

p ■ Z = —m_i_ sinh(?/ — <^), (38) 

where m_L = \/p‘]_ + rn?, y = tanh“^(p^/p°) is the particle’s rapidity, and (p is the particle’s 
azimuthal angle. In order to set up the distribution function, having p^ dehned in Eq. (37), 
one can use (2), (3), (5), and (38) to hnd 

P^E^ipP'' = m± cosh 6± cosh(|/ — — p± sinh6 *_l cos(0 — tp) 

+ ix sinh 6^ cosh(|/ — <^) — cosh 6*j_ cos(0 — p) 

+ sinh^(|/ - g) 

+ iy p\siM{(t)-ip) ■ (39) 

Expanding Eq. (32) for = {t,x,y,z) and contracting with p^, one obtains 


p^d^Tiy = — sin 6d^ 


. / , Md 

p± sm(0 -(p)zri 


A 


dd 


H- m± cos C sin 9 sinh(|/ — <^) ( d cos d + — sin 9 

Ou 

+ cos ((■ sin 6* ( p_i_ cos C, cos(0 — (p) + mj_ sin C, cosh(?/ ~ ‘^) j fd sin ^ ~ ^ ^ 


+ co=c|(px=inCco=(^-<^)-™.co=Cco.h(,-.)) 


d(d(j)d9 , 


( 40 ) 
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where d = d{C,(j),9). 


A. Boost-invariant and cylindrically-symmetric system 


If the system is cylindrically symmetric, one has 

d{C(j),e)^d{ce). (41) 

If, in addition, the system is boost invariant along the beam-line direction, the hypersurface 
in (rr0<^)-space is constant with respect to the longitudinal rapidity <^. As a consequence, 
d{(, 6) sin^ (which is the projection of (i(C, 0) normal to the ?-axis) should be constant, i.e. 

d(C, 0) sin 6 = constant. (42) 


Using this and taking the constant to be the value of the function at a typical longitudinal 
rapidity, i.e. <^ = 0 (0 = 7r/2), one hnds 

d{0 


diC0) = 


sin 9 ’ 


dd{C,e) ^ 


( 43 ) 


89 '"^'^^sin 9 ’ 

where d{() = d{(,9 = 'k/2). Using the above simplihcations, one hnds the boost-invariant 
cylindrically-symmetric form of p^d^'En 


p^d?Y^^ = — (i((C)^ cosC csc^ 6* 


d'{C) sin cos(0 — p) — m± cos ( cosh(|/ — <^) j 


+d{() ( p± cos ( cos(0 — (p) + m± sin ( cosh(|/ — <^) 


dCd(pd9. (44) 


VII. NUMERICAL RESULTS 

We now turn to our numerical results. We present comparisons of results obtained us¬ 
ing the dynamical equations of anisotropic hydrodynamics presented in Sec. IV and the 
second-order viscous hydrodynamics equations from Denicol et ah [29].^ For anisotropic hy¬ 
drodynamics, we use the freeze-out method detailed in Sec. VI and, for second-order viscous 

^ For the smooth initial conditions considered herein, the vorticity is zero at all times. We also set the 
transport coefficient to zero since this has been done in almost all other implementations to date (see, 
however, [35]). As a result, we drop the last two terms in Eq. (D3). 
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hydrodynamics, we use the standard Grad-14 ansatz for the dissipative correction to the one- 
particle distribution function. We present the full details of both the viscous hydrodynamics 
equations and the freeze-out prescription used for our comparisons in App. D. For both 
anisotropic and second-order viscous hydrodynamics we use the lattice-parameterization of 
the EoS presented in Sec. V. Both codes were tested by comparing the evolution of the 
system initialized with a Gubser temperature- and flow-prohle and then comparing with 
the exact solution appropriate to each framework. In both cases, using the lattice spacing, 
temporal step size, etc. specihed below, we were able to reproduce the corresponding exact 
Gubser solution to very high accuracy at all times. We present the details of our code tests 
in App. F. The code used to produce all hgures in this text is publicly available [68]. 

A. Nucleus-nucleus collisions 

For all results presented in this section we use the Glauber wounded-nucleon overlap to 
set the initial energy density. As our test case we consider Pb-Pb collisions with a center of 
mass energy of 2.76 GeV/nucleon. We take the inelastic nucleon-nucleon scattering cross- 
section to be (Ttvat = 62 mb. We use 300 points in the radial direction with a lattice spacing of 
Ar = 0.05 fm and temporal step size of At = 0.01 fm/c. We use fourth-order Runge-Kutta 
integration for the temporal updates and fourth-order centered differences for the evaluation 
of all spatial derivatives.® Unless otherwise indicated, we take the central initial temperature 
to be To = 600 MeV at tq = 0.25 fm/c and assume that the system is initially isotropic, 
i.e. a^iro) = a^(ro) = 1 for anisotropic hydrodynamics and 7r^'^(ro) = 0 for second-order 
viscous hydrodynamics. We take the freeze-out temperature to be = Tpo = 150 MeV in 
all cases. 

1. Hydrodynamic evolution and spectra 

In Fig. 2 we present a comparison of the numerical 1-l-ld solution to the Israel-Stewart 
and anisotropic hydrodynamics equations for the case dvrr^/s = 1. As can be seen from 
Figs. 2a and 2b, for this small value of the shear viscosity to entropy ratio, the two methods 
agree quite well, with only small differences seen in both the temperature and transverse 

® Since the initial conditions considered herein are smooth, naive centered differences generally suffice. 
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FIG. 2. Comparison of the numerical l+ld solution to the Israel-Stewart and anisotropic hydro¬ 
dynamics equations at r = 5.25 fm/c. The shear viscosity to entropy density ratio was taken to 
be 47rr//s = 1. The four panels show (a) the effective temperature, (b) the transverse flow rapidity 
(^±)) (c) the transverse and longitudinal pressures, and (d) the ratio of the LRF longitudinal and 
transverse pressures. 

rapidity profiles at the time shown. At very early times, the differences are somewhat larger, 
but by r = 5.25 fm/c the agreement is quite good. Note that the bumps located at r ~ 10 
are a reflection of evolution near the softest point in the EoS and both the standard viscous 
hydrodynamic method for imposing the EoS and the anisotropic hydrodynamics method 
implemented herein seem to give comparable results. Although we don’t show it here, we 
also performed tests for Airri/s = 0.1 and found that anisotropic hydrodynamics and Israel- 
Stewart second-order hydrodynamics give virtually indistinguishable results in this case. 

In Figs. 2c and 2d we present a comparison of the results obtained for the local rest frame 
(LRF) transverse and longitudinal pressures, again for the case of dvrr^/s = 1. For anisotropic 
hydrodynamics, we compute using Eq. (11) which corresponds to the radial pressure since. 
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(a) 


(b) 



(c) 




PT [GeV] 


PT [GeV] 


FIG. 3. Comparison of the numerical l+ld solution to the Israel-Stewart and anisotropic hydro¬ 
dynamics equations for the same conditions as shown in Fig. 2. The three panels show (a) the 
freeze-out hypersurfaces, (b) resulting neutral pion spectra, and (c) the ratio of the neutral pion 
spectra obtained using Israel-Stewart and anisotropic hydrodynamics. 

using the vectors specified in Eq. (Al), one sees that maps to the radial direction in 
the LRF. In the case of second-order viscous hydrodynamics we identify Pt = Peq + ttJ! 
and Pl = Peq + again in the LRF. As these panels demonstrate, for Airri/s = 1 the 
pressures obtained are similar using both methods with the the maximum difference for 
r < 10 fm being less than approximately 2%. At very large r we see larger differences, 
with the longitudinal pressure predicted by the Israel-Stewart equations becoming negative. 
However, the temperature where this occurs is quite small and far below the freeze-out 
temperature. We do note that one hnds that anisotropic hydrodynamics predicts that the 
system is generally closer to isotropy than Israel-Stewart hydrodynamics. As we will see 
below, the smaller pressure anisotropy effectively reduces the shear stress contribution to 
the transverse pressure, reducing the buildup of radial flow and leading to softer particle 
spectra from anisotropic hydrodynamics. 

Next, we turn to Fig. 3. In this hgure, we present comparisons of the freeze-out hy¬ 
persurfaces (Fig. 3a), the resulting neutral pion spectrum (Fig. 3b), and the ratio of the 
neutral pion spectra obtained using Israel-Stewart and anisotropic hydrodynamics (Fig. 3c). 
As these panels demonstrate one hnds that the freeze-out hypersurfaces and resulting par¬ 
ticle spectra are quite similar for Anrj/s = 1 which should not be overly surprising. As 
Fig. 3c demonstrates the maximal correction for < 3 GeV is approximately 5%, with 
second-order viscous hydrodynamics predicting a harder distribution with a slightly higher 

mean-pr- 
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FIG. 4. Comparison of the numerical l+ld solution to the Israel-Stewart and anisotropic hydro¬ 
dynamics equations. The shear viscosity to entropy density ratio was taken to be A'K'qjs = 3. The 
quantities shown in the panels and the initial conditions are the same as in Fig. 2. 


(a) 



(b) (c) 



Pt [GeV] Pt [GeV] 


FIG. 5. Comparison of the numerical 1-|-Id solution to the Israel-Stewart and anisotropic hydro¬ 
dynamics equations for the same conditions as shown in Fig. 4. The quantities shown in the panels 
are the same as in Fig. 2. 

In Figs. 4 - 8 we present similar plots for the cases Anri/s = 3 and Anri/s = 10. As Fig. 4 


17 












(a) 



(b) 



FIG. 6. Comparison of the vr*^ spectra coming from Israel-Stewart with and without {6f = 0) 
viscous corrections to the distribution function and anisotropic hydrodynamics with a complete 
treatment and assuming isotropic freeze-out using the effective temperature. In this figure, we use 
the same conditions and parameters as shown in Fig. 4. The two panels show (a) the resulting 
neutral pion spectra and (b) the logarithmic slope of the various curves. 

demonstrates, for Anrj/s = 3 we find larger differences between anisotropic hydrodynamics 
and second-order viscons hydrodynamics, as could be expected a priori. One sees that in 
both cases larger pressure anisotropies are generated than for Anrj/s = 1 and the maximum 
difference in the pressure anisotropy approaching 70% in the region r < 10 fm. The freeze- 
out hypersurface and neutral pion spectra shown in Fig. 5 also show larger differences, with 
the final neutral pion spectra having a maximum difference of approximately 18% for < 3 
GeV. Once again, we find that the neutral pion spectra predicted by second-order viscous 
hydrodynamics are harder than that predicted by anisotropic hydrodynamics. 

In order to further understand the differences seen in the predicted neutral pion spectra, 
in Fig. 6 we present the result of using different freeze-out prescriptions for both anisotropic 
hydrodynamics and second-order viscous hydrodynamics. For anisotropic hydrodynamics, 
we show two cases: (1) using the fully anisotropic distribution as in Eq. (39) together with 
the local value of A to construct the local distribution via Eq. (5) and (2) manually setting 

= ^ 1 / = ^2 = 0 in Eq. (39) and setting A —?• T which is the local effective isotropic 
temperature determined from the local energy density. For viscous hydrodynamics, we 
also show two cases: (1) Using the full Grad-14 form given in Eq. (D17) and (2) setting 
6f = 0, which corresponds to discarding the second term in square brackets in Eq. (D17). 
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X [fm] 


FIG. 7. Comparison of the numerical l+ld solution to the Israel-Stewart and anisotropic hydro¬ 
dynamics equations. The shear viscosity to entropy density ratio was taken to be 47rr//s = 10. The 
quantities shown in the panels and the initial conditions are the same as in Fig. 2. 


(a) 



i... 

o- 


Z 

'Tj 


(b) 


(c) 



PT [GeV] 


PT [GeV] 


FIG. 8. Comparison of the numerical l-|-ld solution to the Israel-Stewart and anisotropic hydro¬ 
dynamics equations for the same conditions as shown in Fig. 7. The quantities shown in the panels 
and the initial conditions were the same as in Fig. 3. 

As can be seen from Fig. 6a, for both anisotropic hydrodynamics and second-order viscous 
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hydrodynamics, the inclusion of the (anisotropic) dissipative correction to the one-particle 
distribution function results in a hardening of the particle spectra, however, for anisotropic 
hydrodynamics the correction is slightly smaller. In order to separate the effect of an overall 
shift upwards in the spectra between the various cases, in Fig. 6b we present the logarithmic 
slope of the spectra in all four cases. This plot allows us to more directly see the hardening 
in each case. Even without viscous corrections in the freeze-out distribution function (red 
dashed and green dash-dotted lines in Fig. 6) the spectra obtained using Israel-Stewart 
viscous hydrodynamics are harder than those obtained using anisotropic hydrodynamics. 
This demonstrates that the hardening is primarily due to stronger radial flow, caused by 
larger transverse shear stress and reduced work done by the longitudinal pressure in Israel- 
Stewart hydrodynamics. Viscous corrections at freeze-out further enhance the hardening of 
the spectra obtained using Israel-Stewart hydrodynamics. 

Next we consider the case Anri/s = 10, which is shown in Figs. 7 and 8. As can be 
expected, we see larger differences between the two approaches in this case. At the time 
shown, one sees that the longitudinal pressure predicted by second-order hydrodynamics 
is negative for all r ^ 6.5 fm. If one were to plot the longitudinal pressure predicted 
by second-order viscous hydrodynamics at much earlier times, one would hnd that the 
longitudinal pressure becomes negative in the entire simulation volume. This is indicative 
of the breakdown of second-order viscous hydrodynamics in this case. Not surprisingly, 
as can be seen from Fig. 8, one sees large differences to both the freeze-out hypersurface 
and the hnal neutral pion spectra, with second-order viscous hydrodynamics predicting a 
signihcantly harder distribution and approximately 100% more neutral pions at pt = 3 GeV. 

2. Dissipative particle production 

As our last consideration in the context of Pb-Pb collisions, in Figs. 9 and 10 we present 
the total number of charged particles, Vchg, scaled by the ideal result, Vchg,ideal, as a function 
of dvrp/s. For the purposes of this hgure, we have included the production of charged pions, 
kaons, and protons. In Fig. 9 we took an initial central temperature of Tq = 600 MeV at 
To = 0.25 fm/c and in Fig. 10 we took an initial central temperature of Tq = 500 MeV at the 
same initial time. As both hgures demonstrate, one hnds that second-order viscous hydrody¬ 
namics predicts that A7 h g is a monotonically increasing function oirj/s whereas anisotropic 


20 



FIG. 9. The left panel (a) shows the number of charged particles scaled by the ideal hydrodynamics 
result as function of dvrr^/s obtained using Israel-Stewart viscous hydrodynamics (black line) and 
anisotropic hydrodynamics (red dashed line). The right panel (b) shows the ratio of the Israel- 
Stewart viscous hydrodynamics result to the anisotropic hydrodynamics result. The initial central 
temperature was taken to be Tq = 600 MeV. 



FIG. 10. Same as Fig. 9, except here we take the initial central temperature to be Tq = 500 MeV. 

hydrodynamics predicts that there is a maximum in charged particle production at ATTr]/s ~ 
9-11 depending on the assumed initial temperature. In panel (b) of both hgures we show 
the ratio of the particle production predicted by second-order viscous hydrodynamics and 
that predicted by anisotropic hydrodynamics. For Anrj/s < 5, both Figs. 9 and 10 show 
that the difference in the total number of charged particles produced in Israel-Stewart and 
anisotropic hydrodynamics remains below 10%. However, even a correction on the order 
of 10% could have an important phenomenological impact. Note that the non-monotonic 
behavior as a function of rj/s of particle production in anisotropic hydrodynamics is similar 
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to results obtained in 0+ld anisotropic hydrodynamics [42, 47, 54], which can be straightfor¬ 
wardly understood as the vanishing of dissipative particle production in the free-streaming 

limit. 

B. Proton-Nucleus collisions 

We now consider the case of an asymmetric collision between a proton and a nucleus. Of 
course, since our equations are boost invariant one cannot draw hrm conclusions regarding 
comparisons with experimental results. Our goal is to simply ascertain the magnitude of 
the differences one sees when using anisotropic hydrodynamics versus second-order viscous 
hydrodynamics for small systems. Since the anisotropic hydrodynamics equations have been 
shown to better reproduce the spatiotemporal evolution for small systems subject to Gubser 
flow for all values of p/s [59], we believe that the anisotropic hydrodynamics equations used 
herein should also provide a more faithful reproduction of the bulk evolution for strong 
dynamically generated 1-l-ld flows in small collisional systems at high energies. 

Similar to the previous subsection, we use the Glauber wounded-nucleon overlap to set 
the initial energy density. As our test case, we consider a p-Pb collisions. We use the same 
numerical and physical parameters as in the case of Pb-Pb collisions, except here we consider 
a lower initial central temperature of Tq = 400 MeV. Our hndings are shown in Figs. 11- 
13. From Fig. 11 we see that, assuming 47rp/s = 3, there are relatively small differences 
in the temperature and flow prohles at r = 2.25 fm/c. However, we see quite signihcant 
differences in the transverse and longitudinal pressures, with second-order viscous hydrody¬ 
namics again predicting negative longitudinal pressure and a quite different longitudinal to 
transverse pressure ratio. From Fig. 12 we see that the two methods result in quite different 
particle spectra, even at low momentum. Finally, in Fig. 13, we present a comparison us¬ 
ing different freeze-out prescriptions for both anisotropic hydrodynamics and second-order 
viscous hydrodynamics. The panels and methods are the same as we previously described 
in the context of Fig. 6. As we can see from Fig. 13, the inclusion of viscous (anisotropic) 
corrections to the distribution function has a signihcant effect for p-A collisions. The two 
methods seem to agree qualitatively concerning the direction of the correction, but differ 
quantitatively. 

Gomparing the smaller p-Pb to the larger Pb-Pb collision system, we see that pressure 
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X [fm] X [fm] 


FIG. 11. Comparison of the numerical 1+ld solution to the Israel-Stewart and anisotropic hy¬ 
drodynamics equations at r = 2.25 fm/c. The initial central temperature was Tq = 400 MeV at 
To = 0.25 fm/c with = 3. The quantities shown in the panels are the same as in Fig. 2. 



r [fm] pj [GeV] 


FIG. 12. Comparison of the numerical l-|-ld solution to the Israel-Stewart and anisotropic hy¬ 
drodynamics equations. The parameters are the same as in Fig. 11. The quantities shown in the 
panels are the same as in Fig. 3. 


anisotropy effects are larger in the smaller system, bnt that they are also more severely over- 
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FIG. 13. Comparison of the p-Pb vr^ spectra resulting from Israel-Stewart hydrodynamics with and 
without {6f = 0) viscous corrections to the distribution function, and anisotropic hydrodynamics 
with a complete treatment and assuming isotropic freeze-out using the effective temperature. In 
this figure, we use the same conditions and parameters as shown in Fig. 11. The two panels show 
(a) the resulting neutral pion spectra and (b) the logarithmic slope of the various curves. 

estimated by the Israel-Stewart approach. By resumming the leading viscous effects into 
the leading-order distribution function, as is done in anisotropic hydrodynamics, the pres¬ 
sure anisotropies are signihcantly reduced and the validity of the hydrodynamic approach, 
especially in small collision systems, is signihcantly improved. 

VIII. CONCLUSIONS 

In this paper we demonstrated how to (i) impose a realistic EoS and (ii) self-consistently 
perform hadronic freeze-out in the context of leading-order l-|-ld anisotropic hydrodynamics. 
The methods used here to implement these ingredients can be straightforwardly extended 
to dynamical evolution in more than one spatial dimension. In the case of freeze-out, they 
can be also extended to account for non-ellipsoidal corrections to the LRF one-particle 
distribution function, based on the viscous anisotropic hydrodynamics formalism of Refs. [47, 
53]. 

We compared the results obtained with anisotropic hydrodynamics to results obtained 
with the widely-used Israel-Stewart framework. In the limit of small rj/s we found that 
the two frameworks agree well as they should. For Pb-Pb collisions with Anrj/s = 1, the 
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maximal differences in the pressure anisotropy and neutral pion spectra for pt < S GeV 
were found to be approximately 2%. For larger values of p/s, we naturally saw larger 
differences. For Pb-Pb collisions with Airp/s = 10, we found signihcant corrections to the 
pressure anisotropy and, for pt < 3 GeV, we found a maximal effect of approximately 100% 
for the neutral pion spectra. In general, we found that, for Pb-Pb collisions, second-order 
viscous hydrodynamics predicts stronger radial flow, resulting in harder spectra and more 
dissipative particle production than anisotropic hydrodynamics. 

We also presented results for the total charged particle multiplicity iV^hg as a function p/s. 
We found that, while Whg increases monotonically as a function of p/s using second-order 
viscous hydrodynamics, anisotropic hydrodynamics predicts that there is maximum in iVchg 
around Airp/s ~ 9 — 11 depending on the assumed initial temperature. This observation 
is consistent with the expectation that dissipative particle production should vanish in the 
free-streaming limit. 

Finally, we considered the case of p-A collisions. In this case, we found quite signihcant 
differences between the two frameworks. Since the anisotropic hydrodynamics equations 
used herein have been shown to better reproduce exact solutions to the Boltzmann equation 
in the relaxation time approximation, we have some reason to believe that the anisotropic 
hydrodynamics results obtained herein are more reliable than those obtained using Israel- 
Stewart second-order viscous hydrodynamics. Our results show that viscous effects are 
smaller in anisotropic hydrodynamics than predicted by Israel-Stewart theory, which im¬ 
proves signihcantly the applicability of a hydrodynamic approach to such small collision 
systems. We also demonstrated that the inclusion or exclusion of the viscous (anisotropic) 
corrections to the freeze-out one-particle distribution functions dramatically influences, e.g., 
the hnal neutral pion spectra. This indicates that at freeze-out the system is not equilibrated 
and still quite anisotropic in the LRF. We provided further evidence for this conclusion by 
computing the pressure anisotropy of the fluid for the case of p-A collisions, hnding that 
both the anisotropic hydrodynamics and Israel-Stewart frameworks predict quite large pres¬ 
sure anisotropies even on the late-proper-time portion of the freeze-out hypersurface. In this 
context, we also emphasize that using leading-order anisotropic hydrodynamics one is able 
to guarantee that the pressures and one-particle distribution functions are positive. 

Looking to the future, we have demonstrated how to implement two critically needed 
components for the anisotropic hydrodynamics program. The next steps will be to extend the 
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codes to 2+ld and 3+ld evolution including also perturbatively the second-order corrections 
in the spirit of Refs. [47, 53]. Of course, one could already use the codes developed for use 
in this paper to attempt phenomenological hts of data coming from central Pb-Pb and p-Pb 
collisions at RHIC and LHC. For this, one merely needs to add a hadronic afterburner and 
perform some htting. We postpone this phenomenological exercise to a future publication. 
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Appendix A: Minkowski basis vectors 

Eor azimuthally-symmetric and longitudinally boost-invariant flow, one can parameterize 
the orthogonal basis vectors in lab frame cartesian Minkowski coordinates as 


= cosh 6 *_l cosh , 
= sinh 6*j_ cos 0, 
V? = sinh 6*j_ sin 0, 
= cosh 6 *_l sinh , 


= sinh cosh ^ , 
= cosh 6*_i_ cos 0 , 
= cosh 6^ sin 0 , 
X^ = sinh 0j_ sinh <;■, 


(Al) 



= sinh ? , 

= 0 , 

= 0 , 

= cosh . 


= — sin 0 


COS0 , 


= 0 , 
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In the local rest frame (LRF) they are simply 


■^^LRF ~ (1) 0) 0) 0) ) 

x':j,p = (o,i,o,o), 

^LRF ~ (0) 0) 1) 0) ) 

^(:rf = (0,0,0,1). (A2) 

Appendix B: Formulas for derivatives 

In this section, the identities for derivatives and divergences for 1+ld boost-invariant and 
azimuthally-symmetric flow are summarized. 

Directional derivatives 


Du = u ■ d = cosh 6±dr + sinh 6±dr , 
Dx = X ■ d = sinh 9^dr + cosh 9x_dr , 

D, = Z ■d=-d,. 

T ^ 


Divergences 


9u = d ■ u = cosh6*_L —h dr9±j + sinh0_L -|- dr9±j , 

9x = d ■ X = sinh 9± -|- dr9j^ + cosh 9± , 

5. y = 0, 

9, = d-Z = 0. 


Convective derivatives 


(Bl) 


(B2) 


DuU = {u ■ d)u = X (cosh6*_L(9,-6*_L -|- sinh6'_L(9r6*_L), 

DuX = {u ■ d)X = u (cosh6'_L(9r6'_L -|- sinh6'_L(9r6*_L), 

DuY = (n ■ d)F = 0 , 

DuZ ={u-d)Z = Q. (B3) 
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Directional derivatives 


DxU = (X ■ d)u = X (sinh 6'_i_dT-6'j_ + cosh6'_i_dr6'j_), 
D^X = {X ■ d)X = u (sinh6'_LdT-6'_L + cosh6'_Ldr6'_L), 

= (x-d)y = 0, 

D^Z = {X-d)Z = 


Dyu = {Y ■ d)u = 

DyX = {Y ■ d)X = 

DyY = {Y ■ d)Y = 
DyZ = {Y ■ d)Z = 


sinh^x 

r 

cosh6*j 

r 

- (u sinh 9± 
r 


0 , 


X cosh6*j_), 


n -ry cosh^x 

DzU = [Z ■ o)u = —^— Z , 


DzX = {Z ■ d)X = 


T 

sinh Oj 

T 


DzY = {Z-d)Y = 0, 

DzZ= {Z-d)Z = -(wcosh^j 

r 


X sinh 9 i 


(B4) 


(B6) 


(B6) 


Appendix C: Special Functions 

In this appendix we list the TZ and Ti functions appearing in Sec. IV. They are [48, 59] 


TZ{a^,az) = d(j)a'i'H2 , (Cl) 

T-iTTiax^Ciz) = d(j) al cos^ (j)'H2T ’ (^2) 

V-Tyiax^ctz) = d^OySin^'HaT , (C3) 

1-Ll{oix,oiz) = ^ j d(j)a‘]_'H2L 










with a_L = \ al cos^ (f) + sin^ cj) and 


^2t( 2/) = (^2 _^i)3/2 (|(V - 1) tanh-^ ^^^~^ - y^y^ - 1 j , 

= (^2 ^ 1 ) 3/2 . 


(C5) 

(C6) 

(C7) 


Appendix D: Second-order viscons hydrodynamics 

As with anisotropic hydrodynamics, the viscous hydrodynamics dynamical equations can 
be obtained by taking moments of the Boltzmann equation. Taking the hrst and second 
moments of Boltzmann equation one obtains [10, 29] 

(e + P)D^u^ = V^P - + 7r>^^DuU, , (Dl) 

D^e = -(£ + P)V^u^ + , (D2) 

= 2rja^'^ - , (D3) 

where £ = £eq and P = Peq are the equilibrium (isotropic) energy density and pressure, 
respectively, r^r is the shear relaxation time, and is the shear-shear-coupling transport 
coefficient. The various notations used are 

d^u'' = + r;:„n“, 

Pu = , 

, 

^ (D4) 

The non-vanishing Christoffel symbols for polar Milne coordinates are T^^ = r, T^.^ = 1/r, 
T^^ = —r, and T^^ = 1/r. For the smooth initial conditions considered herein, the vorticity 
is zero at all times. We also set the transport coefficient r^-Tr to zero since this has been done 
in almost all other implementations to date (see, however, [35]). As a result, we drop the 
last two terms in Eq. (D3). 
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1+ld equations of motion 


In the boost-invariant and cylindrically-symmetric case, one has = {u^, m'’, 0, 0) and, as 
a result, v = tanh6'j_ = ju^. In addition, for this case, the shear tensor has the following 
form 


/0 0 \ 

q q 

0 0 TT'^^ 0 

0 0 0 


(D5) 


Note also that, in the boost-invariant and cylindrically-symmetric case, Du as dehned in 
Eq. (D4) reduces to the expression presented previously in Eq. (Bl). 

In this case, expanding Eqs. (Dl), (D2), and (D3) in polar Milne coordinates one obtains 
Eve independent equations 


{e + P)DuU^ = {drP - d.O - {drP - d ^), (D6) 

{e + P)DuU^ = -u-^vT {drP - dX) - {drP - d^O , (D7) 

Due = -{e + P)eu - KY - (D8) 


where 


= -2ri - ttJ , 

T^{Du7il + ‘^Ou'nl) = -2r] vX - ttJ , 


- di,7r^ = V drill + vdrlll nl drV +drV-i -H , 


V V 


T r J T 

dylll = V drill + drill + nl ( drV + - + - --TT^ . 


r 


In addition, one needs the following 


_ u'^DuU'^ + -{uY^Ou , 

o 

r2v<X = -- + 

r 6 

rXX = -- + id„, 

r 3 

?y^ ?7^ 

9u = 'VaU‘^ = daU^ = drVd' -|- drVd' H-1-, 

r r 


(D9) 

(DIO) 

(Dll) 

(D12) 

(D13) 

(D14) 

(D15) 

(D16) 
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where = —vti^ and = —vr^ ~ (1 ~ which are a consequence of the transversality 
of the shear-stress tensor, = 0. This system of equations has to be closed by providing 

an equation of state (EoS), e.g. Peq = Peqi^eq)- For the numerical results presented in the 
body of the manuscript, we use the lattice-based EoS specihed in Sec. V. 


Viscous hydrodynamics freeze-out 

The distribution function on the freeze-out hypersurface can be computed assuming that 
there is a linear correction to the equilibrium distribution function [69] 


f{p,x) = /, 


eq 


1 -|- (1 — a/eq) 




flip 


2{e + P)T^ 


(D17) 


Using tensor transformations applied to Eq. (37) the components of the four-momentum in 
polar Milne coordinates are 


= p* cosh ^ — p^ sinh = m± cosh(p — <^), 
p^ = cos (j) + p'^ sin <p = p± cos(0 — ip ), 
p^ = —p-h p^-=-sm(0 — ip ), 

f 

^ ,sinh<^ .cosh<^ mj_ . , , , 

p^ = —p-h p - = -smh(p — . 

T T T 


Using Eq. (D5) and expanding p^Pi/Vr^^ in polar Milne coordinates one has 


(D18) 


PiiP^TT'^ = Pp” + PtPtTt" + PrPrTT" + p;!r" + p^n** + p\ 


TT . 


(D19) 


As a result, one obtains 


=- 1 

r0^2 „• 2 


m±v cosh(?/ — — p± cos(0 — <p) j 


TT^p^ sin^(0 -ip) - vrjm^ sinh^(p - <^) . 


(D20) 


Appendix E: Factorization in non-conformal anisotropic hydrodynamics 

In the case that the particles comprising the system are massive, conformality is broken 
and it is no longer possible to multiplicatively factorize the energy density and pressures as 
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m mm 


FIG. 14. The ratios (a) r^, (b) tt, and (c) ri as a function of m = m/X for three different values 
of the pressure anisotropy Pl/Pt = {0.1,1,10} corresponding to = {0.65499,1, 5.4631}. 


in Eqs. (10)-(13) [42, 49, 51, 56]. For a system with constant mass m one has instead [51] 

£ = 7 / 3 ( 0 :, <h, rh) , (El) 


Pt = 'H3T{oL,^,m) X'^, 
PL = 7/3L(a,$,m) A^ 


(E2) 

(E3) 


where o = ay, a^), m = m/X, p = p/X, and the constraint 4> = ^ ~ 1 with 

i G {x, y, z} is implicit. The functions 7 / 3 , PLst, and PL^l are given by 

'H3{cx,^,m) = Na^ay [ d(l)a]_ f 7/2 (^ 4 ) 


- 27 r 


Q:_l a:_Lp / 

rh 


'H 3 T{oL,^,m) =-Na^ay d(j)a]_ dpp^feJ\/p^ + m^') PL 2 t{—, ) , (E5) 

Jo Jo ^ ^ \Q:± Oi±p^ 


1*277 1*00 . \ ( OL 

'H 3 Lict,^,m)=Naxay d(j)a]_ dpp^feq (^\/ p^ + dpj 7/2Lf ^, 


m 


\q:_l a:_Lp 


(E6) 


with cos^ 0 + a^sin^^, p = |p|, N = Aldof/(27r)^ with being the number of 

degrees of freedom, and the functions 7 / 2 , JJ-2T, and 7 / 2 L are given by 


P2{y,z) = 


Ji2T{y,z) = 


P2L{y,z) = 


y 


Vy 


2 - 1 
y 


(1 + z^) tanh 


-1 y 


2 _ 1 


1/2 + ^2 


+ V(y'^ + ^^)iy‘^ - 1) I > 


(//2 - 1 ) 3/2 

y^ 

(//2 - 1 ) 3/2 


( 2 ;^ + 2y‘^ — 1 ) tanh ^ 


2/2-1 


2/2 + 2;2 

\/ (2/2 — 1)(2/^ + 2^2) — (2;2 -g 1) tanh“^ 


- Viy‘^ - i)(i/^ + ^^) 


y^-i 
2/2 + ^2 


(E7) 


, (E8) 


(E9) 


Since 7 / 3 , PLot, and PLsl depend on many variables, we first restrict ourselves to the case 
that Qfa; = tty which is appropriate for the case of 1+ld dynamics considered in the body of 
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the text. In addition, since, in the body of the text we assume $ = 0, we also assume this 
here. With these assumptions, one can write in terms of [51] 

2a? 


a^ =-^ 

" Sa? - 1 


(ElO) 


In order to test the degree to which the conformal factorization occurs in the massive case, 
we dehne three ratios 

'H3(Q:,0,m) 


r. = 


Tt = 


ri = 


'H3T{l,0,rh)'HTx{oix,Oiz) ’ 
'H3L{cx,0,m) 


(Ell) 

(E12) 

(E13) 


'H3L(l,0,m)'HL(aa;,a^) ’ 
where TZ, T-Ltx, and "Hl are dehned in Eqs. (C1)-(C4). The factors of 'H 3 (l, 0, m), 'H 3 'r(l, 0, m), 

'HL(l,0,m) are introduced in the denominator in order to compensate for the trivial mass 
dependence of the EoS in the isotropic case. If these ratios are equal to one, then there is 
perfect factorization of the diagonal components of the energy momentum tensor and the 
size of the deviation from unity is indicative of the degree to which factorization is broken 
in the non-conformal case. We present our numerical evaluation of these three ratios in 
Fig. 14. As can be seen from the three panels, as long as 0.1 < Pl/Pt < 10 and m < 1, 
the maximum correction to all of these ratios is approximately 5%. For the case of oblate 
pressure anisotropies, which is more relevant for phenomenological application to heavy ion 
collisions, the maximum correction is approximately 3%. 


Appendix F: Gubser flow tests 

In this appendix we present results of code tests for both the anisotropic hydrodynamics 
and viscous hydrodynamics codes used in the body of the text. For our tests, we initialize the 
codes at a given proper time using the exact Gubser solution of the hydrodynamic equations 
appropriate for each case (see [59, 70]).® Gubser flow is a conformal flow which, in polar 
® These exact solutions are the ones appropriate for each approximation scheme and are not to be confused 
with e.g. the exact solution of the Boltzmann equation in the relaxation time approximation (RTA) 
[57, 58]. When compared to the exact RTA solution, anisotropic hydrodynamics does a dramatically 
better job reproducing the exact RTA solution and can be shown to analytically reproduce both the ideal 
and free-streaming limits [59]. 
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FIG. 15. Comparison of the numerical l+ld solution to the second-order viscous hydrodynamics 
equations with the corresponding exact Gubser-flow solution [70] at r = 10 fm/c. The initial 
central temperature was Tq = 600 MeV at tq = 1 fm/c with Air-q/s = 3. For this test we took 
g = (1 fm/c)“^. 

Milne coordinates, is completely determined by symmetry constraints to be [71, 72] 

= (cosh 0^, sinh6'_L, 0, 0), (FI) 


where 


= tanh 


-1 


2q^Tr 


(F2) 


1 + + g2^2 

with q being an energy scale which sets the transverse spatial size of the system. Gubser 
flow is best understood through the introduction of de Sitter variables p and 6 [71] 

1 - + g^r^ 


sinh p = 
tan^ = 


2gr 

2gr 


(F3) 

(F4) 


1 + g^r^ — g^r^ ’ 

where p and 6 are components of de Sitter coordinates, = (p, 6^, 0, <^). Note that, for hxed 
r, the limit r —)■ O’*" corresponds to the limit p —)■ —cx) and the limit r —>■ oo corresponds to 
the limit p —)■ cxd. As a consequence, the de Sitter map —oo < p < +cxd covers the future 
(forward) light cone. In what follows in this appendix, all Weyl-rescaled variables dehned 
in de Sitter coordinates are indicated with a hat. Note that, since Gubser flow is conformal, 
the EoS used in the tests presented below is a conformal (ideal) equation state. 


Gubser flow using Israel-Stewart second-order viscous hydrodynamics 

For the second-order hydrodynamic approximation subject to Gubser flow one has to 
solve two coupled ordinary differential equations subject to a boundary condition at p = po. 
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For the case of Israel-Stewart viscous hydrodynamics, the necessary equations are [70] 

1 dT 2 1 - c/ \ 1 

3 tanhp = -7rJ(p) tanhp, (F5) 

4 o 4 

+ - (tt^) tanh p + ^ = — tanh p , (F6) 

CL p o '^TT -L 0 

where = iryiTs) and = 5p/{sT). Once the solution of these ordinary differential 
equations is obtained, one can map them back to Minkowski space using Eqs. (F3) and 
(F4). 

We present our results for the comparison of our Israel-Stewart solver and the exact 
solution in Fig. 15. We use the same algorithm, lattice spacing, and temporal time step as 
in the main body of the text. As one can see from Fig. 15, we are able to obtain excellent 
agreement between our numerical solution of the 1-l-ld Israel-Stewart partial differential 
equations and the exact solution subject to Gubser flow, even as late as r = 10 fm/c. 
There are some small discrepancies near the boundary of the simulated region which are 
due to boundary effects. We have checked that these effects can be reduced by using larger 
simulation volumes. 


Gubser flow using anisotropic hydrodynamics 


The dynamical equations needed to describe the de Sitter-space evolution of a system 
subject to Gubser flow using anisotropic hydrodynamics are [59] 


+ tanh p f ^ + 2 I = 0 , 


dp 


3d2 - 1 


dp 


n^iv) 


6a, da, 3 (3d^ - + i) (f\ ^ ^ 

( X I + tan p - , 


(F7) 

(F8) 


where y = a,lae = ^/{^af--V)J2. The "H-functions appearing above are dehned in 
Eqs. (C5)-(C7). The set of equations can be closed by using the dynamical Landau matching 
condition 




1/4 


A. 


(F9) 


Once the solution of these ordinary differential equations is obtained, one can map them 
back to Minkowski space using Eqs. (F3) and (F4). 
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FIG. 16. Comparison of the numerical l+ld solution to the anisotropic hydrodynamics equations 
with the corresponding exact Gubser-flow solution [59] at r = 10 fm/c. The parameters, grid 
spacings, etc. used are the same as in Fig. 15. 

We present our results for the comparison of our anisotropic hydrodynamics solver and 
the exact solution in Fig. 16. As before, we use the same algorithm, lattice spacing, and 
temporal time step as in the main body of the text. As one can see from Fig. 16, we are 
able to obtain excellent agreement between our numerical solution of the l+ld anisotropic 
hydrodynamics equations and the exact solution subject to Gubser flow, even as late as 
r = 10 fm/c. As in the case of the Israel-Stewart solver, there are some small discrepancies 
near the boundary of the simulated region which are due to boundary effects. We have, once 
again, checked that these effects can be reduced by using larger simulation volumes. 
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